#!/usr/bin/perl
# 2013/09/16 distance vs association frequency calculation to check self ligation

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [map file] -o [output file]\n";
}

#50:AGGTC        II      2142721 -       37      201160R 2142719 I       125609  -       37      100065L 125603  U:U     1:1     NA:NA
#51:ATTTC        II      2160888 -       37      201172R 2160886 II      2162163 -       37      201174L 2162157 U:U     1:1     NA:NA
#52:GGGGG        II      3400768 -       37      201873R 3400766 I       4177332 -       37      102398L 4177325 U:U     1:1     NA:NA
#55:TATTT        I       3677215 -       37      102141R 3677213 I       1615581 +       37      100966R 1615579 U:U     1:1     NA:NA
#56:GAGGG        II      3718799 -       37      202032R 3718797 II      3718803 -       37      202032L 3718797 U:U     1:1     NA:NA
#57:TACAA        I       2592440 -       37      101553R 2592438 II      3315132 +       37      201825R 3315130 U:U     1:1     NA:NA


my %opt;
getopts("i:o:", \%opt);
my $FILE_map = $opt{i};
my $FILE_out = $opt{o};

my $Resolution = 100;

#---------------------------------------
# collect information
#---------------------------------------
my %data;
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;

	if($chr1 ne $chr2){
		next;
	}

	if($type ne 'U:U'){
		next;
	}

	my $hinNum1 = substr($hinID1, 0, -1);
	my $hinNum2 = substr($hinID2, 0, -1);
	my $hinSide1 = substr($hinID1, -1);
	my $hinSide2 = substr($hinID2, -1);

	my $distance= abs($hinLoc1 - $hinLoc2);
	my $category = int($distance / $Resolution) * $Resolution;

	# definition of association
	# R---1---L R---2---L
	# 1R - 2L   out-out
	# 1L - 2R   in-in
	# 1R - 2R   out - in
	# 1L - 2L	in - out

	my $AsociationType = '';
	if($hinNum1 < $hinNum2){
		if($hinSide1 eq 'R'){
			if($hinSide2 eq 'R'){
				$AsociationType = 'out-in';
			}elsif($hinSide2 eq 'L'){
				$AsociationType = 'out-out';
			}
		}elsif($hinSide1 eq 'L'){
			if($hinSide2 eq 'R'){
				$AsociationType = 'in-in';
			}elsif($hinSide2 eq 'L'){
				$AsociationType = 'in-out';
			}
		}
	}else{
		if($hinSide2 eq 'R'){
			if($hinSide1 eq 'R'){
				$AsociationType = 'out-in';
			}elsif($hinSide1 eq 'L'){
				$AsociationType = 'out-out';
			}
		}elsif($hinSide2 eq 'L'){
			if($hinSide1 eq 'R'){
				$AsociationType = 'in-in';
			}elsif($hinSide1 eq 'L'){
				$AsociationType = 'in-out';
			}
		}
	}

	if($AsociationType eq ''){
		die "Association type is not defined for the data\n$_\n";
	}

	$data{$category}{$AsociationType}++;
}
$fh_map->close();


#---------------------------------------
# output distance vs freqeuncy
#---------------------------------------
my $fh_out = IO::File->new($FILE_out, 'w') or die "cannot write $FILE_out: $!";
foreach my $distance(sort{$a <=> $b} keys %data){
	$fh_out->print("$distance");
	foreach my $type(qw(out-out in-in out-in in-out)){
		if(exists $data{$distance}{$type}){
			$fh_out->print("\t$data{$distance}{$type}");
		}else{
			$fh_out->print("\t0");
		}
	}
	$fh_out->print("\n");
}
$fh_out->close();